#include <math.h>

#include <vector>

#include "Interpolation2D.h"
#include "Interpolation3D.h"
#include "function.h"

void test1() {
  int m = 7, n = 7;
  float xa[] = {0.00117378, 0.00472397, 0.0171859, 0.0409838,
                0.0780596,  0.130081,   0.198535};
  float ya[] = {0.365616,  1.893040,  3.786090, 7.572170,
                15.144300, 30.288700, 60.577400};
  float za[7][7] = {0.0217822, 0.0253224, 0.0288237, 0.0346827, 0.0448323,
                    0.0636086, 0.100366,  0.0233179, 0.0268545, 0.0303556,
                    0.0362159, 0.0463659, 0.0651426, 0.101902,  0.0296429,
                    0.0331470, 0.0366371, 0.0425000, 0.0526603, 0.0714467,
                    0.108208,  0.0402311, 0.0440292, 0.0477457, 0.0538394,
                    0.0641187, 0.0829203, 0.119654,  0.0511250, 0.0554077,
                    0.0595859, 0.0662932, 0.0771901, 0.0963434, 0.133061,
                    0.0625876, 0.0673198, 0.0719785, 0.0794046, 0.0910973,
                    0.110757,  0.147656,  0.0748282, 0.0800098, 0.0851434,
                    0.0933663, 0.106111,  0.126669,  0.163872};
  float *zp[7];
  for (int i = 0; i < 7; i++) zp[i] = &za[i][0];
  SplInt2D<double, float> p = SplInt2D<double, float>(m, n);
  p.construct(xa, ya, zp);
  std::vector<float> xv = {0.0001, 0.0008, 0.006, 0.02, 0.12, 0.25};
  std::vector<float> yv = {0.1, 0.5, 2.0, 10.0, 45.0, 100.0};
  for (float y : yv)
    for (float x : xv)
      printf("x = %.4f, y = %.4f, f(x, y) = %f \n", x, y, p.evaluate(x, y));
}

void test2() {
  int m = 3, n = 3, l = 3;
  float *xa = new float[m];
  float *ya = new float[n];
  float *za = new float[l];
  float ***va = new float **[m];
  for (int i = 0; i < m; i++) xa[i] = i;
  for (int j = 0; j < n; j++) ya[j] = j;
  for (int k = 0; k < l; k++) za[k] = k;
  for (int i = 0; i < m; i++) {
    va[i] = new float *[n];
    for (int j = 0; j < n; j++) {
      va[i][j] = new float[l];
      for (int k = 0; k < l; k++) {
        va[i][j][k] = plinear(xa[i], ya[j], za[k]);
        printf("x = %d, y = %d, z = %d, va(x, y, z) = %.2f \n", i, j, k, va[i][j][k]);
      }
    }
  }
  QuaInt3D<double, float> p = QuaInt3D<double, float>(m, n, l);
  p.construct(xa, ya, za, va);
  p.printCoefficients();
  for (int i = 0; i < m; i++) {
    for (int j = 0; j < n; j++) {
      for (int k = 0; k < l; k++) {
        float x = (float)i;
        float y = (float)j;
        float z = (float)k;
        if (p.evaluate(x, y, z) != va[i][j][k])
        {
          printf("x = %d, y = %d, z = %d, f(x, y, z) = %.2f \n", i, j, k, p.evaluate(x, y, z));
          printf("x = %d, y = %d, z = %d, va(x, y, z) = %.2f \n", i, j, k, va[i][j][k]);
        }
      }
    }
  }


}

int main() { 
    test2(); 
}
